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We present a new hydro code based on spectral methods using spherical coordinates. The first 
version of this code aims at studying time evolution of inertial modes in slowly rotating neutron 
stars. In this article, we introduce the anelastic approximation, developed in atmospheric physics, 
using the mass conservation equation to discard acoustic waves. We describe our algorithms and 
some tests of the linear version of the code, and also some preliminary linear results. We show, in 
the Newtonian framework with differentially rotating background, as in the relativistic case with 
the strong Cowling approximation, that the main part of the velocity quickly concentrates near the 
equator of the star. Thus, our time evolution approach gives results analogous to those obtained by 
Karino et al. [Mwithin a calculation of eigenvectors. Furthermore, in agreement with the work of 
", Lockitch et al. J2j, we found that the velocity seems to always get a non- vanishing polar part. 
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In 1998, Andersson || discovered that in any relativistic rotating perfect fluid, r-modes 1 are unstable via the 
coupling with gravitational radiation. This is just a particular case of the Chandrasekhar-Friedman-Schutz (CFS) 
mechanism discovered in 1970 (see j(| and Q). However, a characteristic of r-modes is the fact that this instability 
\ is expected whatever the speed of rotation of the star (see also Friedman et al. 1998 ||). This discovery triggered 
■ an important activity in the neutron stars (NS) community since such inertial instabilities might contribute to 
\Q | explanation of the observed relative slow rotation rates of young and recycled NS (see Bildsten Q). Moreover, such 
an instability could in principle make of NS efficient sources of gravitational waves (GW) for observation by the 
ground interferometric detectors under construction or their direct descendants. Nevertheless, as various not always 
well-known physical processes take place which might inhibit the growth of r-modes in real NS, their relevance is 
still questionable at this time. More details can be found in the reviews by Andersson and Kokkotas 2001 [jl0| and 
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Friedman and Lockitch 2001 |1 

The present study was originally motivated by the results announced by Lindblom et al. in 2001 E2| . These 
authors computed, in the highly non-linear regime, the evolutionary track of the CFS instability of r-modes in a 
^Jy Newtonian NS with a magnified approximate post-Newtonian radiation reaction (RR) force. In their calculation, 
they found a rapid growth of the mode, until strong shocks appeared and quickly damped it. A recent article by 
Arras et al. |ll| asserts that this phenomenon was just an artifact linked with the huge RR force. Yet, we found 
those results so interesting and so surprising that we decided to try to reproduce and better understand them with a 
different approach. The code of Lindblom et al. ]l2[ was written in cylindrical coordinates with 3D finite difference 
scheme and they studied fast rotating stars. We chose to create a spectral code in spherical coordinates and to begin 
with a slowly rotating NS. 

The use of spectral methods was motivated by the fact that they can provide a deep insight in describing the 
turbulence generated by quadratic terms. This is analogous to what is done with Fourier analysis in the study of 
homogeneous turbulence (e.g. Lesieur 1987 Furthermore, we chose a slowly rotating NS for several reasons 

that will be explained in the following. 

Since NS are known for their quite rapid rotation and since inertial modes have for restoring force the Coriolis 
force, an a priori reasoning would probably lead to the conclusion that the only - or most - interesting case of 
inertial modes to study is the case of modes in a fast rotating NS. However, the final answer is not so easy to decide, 
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1 In this paper, we shall use the terminology of the hydro community : inertial modes are the oscillatory modes of a fluid whose restoring 
force is the Coriolis force. .R-modes [or purely toroidal (axial) modes] then form a sub-class of inertial modes. The later were discovered 
by Thompson in 1880 H] . The reader can find a short point of history in Ricutord 2001 H . 
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as among observational data and among models of cooling of pulsars, many things in fact support slowly rotating 
newborn single NS. For instance, the estimations of the rotation periods at birth of the two historical X-rays and 
radio pulsars whose ages are known exactly (the 820 year old 65.8 ms period pulsar in SNR 3C58 Jl6| and 

the 948 year old 33 ms period Crab pulsar) are respectively 60 ms [JL5| and 14 ms |l7|. Moreover, compared with 
assessment of the angular momentum of isolated supergiant stars (whose core collapses will give type II supernovae 
and then potentially NS), these numbers show that important losses of angular momentum are expected to happen 
during the stellar evolution. Furthermore, this is the same concerning the giant phase of less massive stars that give 
white dwarfs. Indeed, the observed rotation periods of 19 white dwarfs range between 12.1 min and 12 days with a 



median value of about 1 hour |18 but, if the sun shrank without losing any angular momentum into a white dwarf of 
the same mass with a radius of 5000 km, its rotation period would be of a few minutes only. In the stellar evolution 
community, several ways to explain these weak angular momenta can be found and the final answer is still not clear. 
But a quite admitted idea is that magnetic field plays probably the key-role via magnetic braking type mechanisms 
(Schatzman 1962 [fl9|). And whatever the actual mechanism, the current conclusion in this community is that a 
huge amount of angular momentum is supposed to be lost during the stellar evolution for main sequence stars of any 
mass. Then, the most difficult thing to explain in stellar evolution physics seems to be the reason why these rotation 
periods at birth (of the orders of 1 hour for white dwarfs and of 1 ms for NS) are so small (e.g. Spruit 1998 |2(| and 
Spruit and Phinney 1998 J2l]]): stellar evolution scenarios do not expect of baby NS to be fast rotating. 

However, it can be mentioned that some isolated fast rotating NS are found. But, for the most part of them, 
the current idea is that these stars have been recycled in binary systems. In these systems, where fast rotating 
NS exist, the NS is thought to have been spun up by the accreted matter. The recent discoveries of accreting 
millisecond pulsars (XTEJ0929-314, SAXJ1808.4-3685, XTEJ1751-305 and 4U 1636-53) whose rotation frequencies 
are respectively 185 Hz, 401 Hz, 435 Hz and 581 Hz (corresponding to respective periods of 5.4 ms, 2.5 ms, 2.3 ms 
and 1.7 ms) support the above scenario. Yet, there is also an example of single ms pulsar associated to a supernova 
remnant: PSR J0537-6910 whose period is 16 ms, whose age is about 5000 yr and which is supposed to have had an 
initial period of a few ms (see for instance [^3|). Hence, in spite all the above arguments, the existence of rapidly 
rotating baby NS should not be completely excluded. Moreover, a fast rotating baby NS with very weak magnetic 
field and consequently very weak losses of angular momentum due to magnetic braking like mechanisms during the 
supergiant phase would escape to the observations as pulsars but would be very interesting as GW sources. 

Nevertheless, such a discussion does not make clear the important point. Indeed, whatever the value of the 
frequency, it does not directly tell if a pulsar is fast rotating or if it is not. What says, in a particular work, if a NS 
can be regarded as (quite) slowly rotating is the relative importance of the deformation for this study. In a more 
general framework, this question is settled by the ratio between the angular velocity of the pulsar and its Keplerian 
angular velocity. Even for the fastest rotating known pulsars (that are in binary systems) this ratio is less than a 
third, since the Kepler frequency is around 1 ms (sec |Q). Furthermore, what is implied in the hydrodynamics of a 
NS is not this ratio, but its square. Hence, in appropriate units, this factor is less than ten per cent of the Coriolis 
force for most of the NS. Anyway, there is a last crucial issue. In a Newtonian star, even if this factor is a few per cent 
correction in the equation of motion, it is fundamental since it creates a coupling between polar and axial part of the 
velocity. Yet, in a relativistic star, the situation is completely different as noticed Kojima psj . Indeed, in a relativistic 
rotating star, the frame dragging term has the same qualitative result: it makes a coupling between the polar and 
the axial parts of the velocity. But, in appropriate units, this coupling is scaled by the ratio between the angular 
velocity and the Keplerian angular velocity and not by the square of this ratio. Hence, even for the fastest rotat- 
ing known NS, the deformation introduces a kind of second order correction that can be neglected in a first approach 2 . 

Thus, our choice of using the slow rotation limit as a first step was motivated by all the astrophysical reasons 
mentioned above: even for ms pulsar, the slow rotation approximation is still quite good. But, secondly, we wanted to 
look for a possible saturation due to non-linear coupling that may occur before a highly non-linear regime is reached. 
To better understand such a phenomenon, we thought it was probably wiser to begin with an easier situation in 
which there are not several effects with consequences of the same orders of magnitude. Thus, we began to build a 
non-linear hydrodynamics code using the Newtonian theory of gravity and the slow rotation approximation. But, 
once a linear Newtonian code had been written (first step to a non-linear version), upgrading it to a general relativity 



2 As an example, we verified that, for a rotation frequency of 300 Hz in a star of 1.7 Mq with a fully relativistic code for stationary 
configuration of rotating stars, the coupling between the spherical harmonic terms due to the drag effect is one order of magnitude larger 
than the coupling due to the deformation of the star. 
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(GR) linear code with strong Cowling approximation 3 was quite obvious and we chose it for a second approach to 
inertial instabilities. Indeed, as it was already mentioned above, Kojima first noticed in 1998 that the frame 
dragging phenomenon makes the relativistic r-modes quite different from the Newtonian r-modes. But concerning 
this linear relativistic study, we also decided to begin with the slow rotation limit to try to get a better understanding 
of the spherical relativistic case, before putting what can be seen as a second order correction. 

This article is organized as follows: in Section ||, we start by describing the physical conditions we chose in the 
Newtonian case. Then we write the Navier Stokes equations (NSE) in dimensionless form, and the RR force we 
introduced in the linear study. This enables us to calculate the associated characteristic numbers and to observe 
that the corresponding numerical problem is a stiff one with typical time scales being orders of magnitude different. 



We discuss some approximations that can be made in Section III. The most important of them is the anelastic 



approximation. Section IV explains the basic tests of the code, with for instance the linear I = m r-modes in a 
rigidly rotating spherical fluid. This admits an easy and analytical solution of Euler equation (EE) in the divergence 
free approximation. We demonstrate that, thanks to spectral methods, the velocity profile is exactly determined and 
preserved within the round-off errors. We also show that the error in the conservation of the energy is only due to 
the time discretization and that it vanishes as the cube of the time-step At. Finally, this section ends with tests of 
the anelastic approximation and of the RR force we adopted. Section |y| deals with results obtained in a differentially 
rotating NS. As expected, pure r-modes are replaced with inertial modes which are still CFS unstable, even if they 
are no longer purely axial. Moreover, we show that if the amplitude of the difference of angular velocity between 
the interior and the surface of the star is sufficiently large, these modes quite rapidly develop a kind of "singularity" 
in the first radial derivative of the velocity. There is then a kind of concentration of the motion near the surface 
and mainly near the equatorial plane. Our time evolution approach thus provides results of the same kind of those 
obtained by Karino et al. in 2001 JjJ with a modes calculation. A viscous term has to be added in order to regularize 
the solution, but it does not change the qualitative result. We shall discuss this result in the conclusion. In Section 



VI, we give some first results for the GR case in a slowly rotating NS with the strong Cowling approximation. In this 



framework, the above conclusions still hold. More results in GR will be given in a following article in preparation. 



The Section VII contains the conclusion and discussion. A detailed description of the numerical algorithm is given in 
the appendix. 



II. EQUATIONS AND NUMBERS 
A. The barotropic case in Newtonian gravity 

Cooling calculations (e.g. Nomoto and Tsuruta 1987 j2Sf| , Yakovlev et al. 2001 J3(J) showed that several minutes 
are enough to enable the NS matter to fall far below its Fermi temperature, that is roughly 10 10 K. Furthermore, in 
the Newtonian non- linear hydrodynamics of a not too young slowly rotating NS, it is not worth trying to use a very 
sophisticated EOS. We shall therefore adopt a barotropic EOS. With those assumptions, we have 

P = P[n], (1) 

where P is the pressure and n the mass density, and the Poisson equation for the gravitational potential U is 

AU = AnGn. (2) 

The Newtonian Navier-Stokes equations (written in the inertial frame for a rigidly rotating NS) are 

(d t + nd v ) W+ (W-W)W+ 2QAW + VH= (3) 
i (v(Vr/-VF) +VA (wAVrj) + r/AW - WAn + V {(( + §) V . w"j\ + A. 

Here we take (r, i9, (p) for a spherical system of coordinates in the inertial frame. In this system, ■& — coincides 
with the direction of the rotation axis of the NS, which is parallel to the angular velocity: Vt = ile z . W is the velocity 



3 In December 2001, during a workshop on r-modes which took place at the Meudon site of Paris Observatory, Carter suggested to call 
strong Cowling approximation the approximation in which all the coefficients of the metric are frozen. This name was chosen to contrast 
with what should be called the weak Cowling approximations where some perturbations of the metric are allowed. See for instance the 
work of Ruoff et al. pa], J27I and Lockitch et al. B using the Kojima equations |Eq|. 
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in the corotating frame, i.e. the part that is added to the velocity of the rigidly rotating background when a mode is 
present. Note that both velocities can be of the same order in the non-linear case. This is the reason why we shall 
not use the term Eulerian perturbation for W . Otherwise, H = J ^^"dn is the enthalpy, rj and £ are respectively 

the dynamical shear and bulk viscosity coefficients, and finally A contains any external accelerations (or force per 
unit of mass) . The modifications that are needed in the case of the linear study with differential rotation or in the 
relativistic linear case are respectively discussed in sections [v| and [vj. 

In what follows, A is mainly the effective gravitational acceleration, i.e. the gradient of the difference between the 
centrifugal potential \£P p 2 and the gravitational potential U, where p is the distance from the rotation axis. In the 
linear regime, we will sometimes introduce a RR acceleration (cf. Blanchet 1993 and 1997 pl|) of the form 
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with (km) =i(fcm + mfc), v l being the full velocity and the superscript ^ the fifth time derivative, which can not be 
easily calculated in a numerical work (see Rezzolla et al. 1999 p2[). Note that this formula is valid only if written 
with Cartesian components of the tensors and this is the reason why we did not really make distinction between 
contra and covariant components. Finally, we insist on the fact that such an acceleration does not include any mass 
multipole but only the current quadrupole that is the most important coefficient for the emission of GW by axial 
modes in a slowly rotating NS. 



To fulfill the system of equations, we need to add the mass continuity equation: 

(d t + nd v ) n + W- Vn + nV- W = 0. 



(5) 



Taking into account that the EOS is a barotropic one, the preceding equation can be written in a slightly different 



way, but we will discuss it in the Section III 



As far as inertial modes are concerned, the typical time scale and length scale are 27rO _1 and R, characteristic 
length of the background star. The velocity associated with those values, R f2, can be very far from the characteristic 
velocity of the mode. Another velocity, scaling that of the mode, must be introduced. But instead, we introduce a, 
a pure number that is defined as the ratio of these two velocities. Otherwise, the typical mass is obviously the star's, 
M* ~ 1.4M©. Bearing all this in mind, we define the following dimensionless variables: 
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Thus the motion equations are written [with same conventions than in Eq. (^)] as 

(d T + d v ) V + a (V ■ V) V + 2e z A V + Vff = 
| (v (Vs • 9) + V A (V A Vs) + sAV - VA s + V ^(A b + §) V 
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where the V and A operators are now performed with dimensionless variables (£, <p) instead of (r, 

Here, we have introduced the following notations: 

- a a pure number, characteristic of the initial amplitude of the mode, but also of the ratio of the non-linear term 
and the Coriolis force. With our conventions, a is twice the usual Rossby number and V is a vector whose norm 
at t — is equal to 1 in a given point on the equator; 

- rio = M+/(^irR 3 ). For the full slow rotation limit approximation (only terms linear in fl and spherical shape), 
R would be the radius of the star and then n the mean density. 

- s(hear) and 6(ulk) dimensionless functions, and v and A pure numbers, all chosen in such a way that if there 
is any viscosity, s and b are equal to 1 in a specific position. Typically, for a single fluid model, this would be 
the center of the star. Nevertheless, it is worth pointing out that to build a more realistic model of NS, several 
different layers could be made to coexist. In this case, there would be different v and A, hugely depending on 
the shell (see for instance Haensel et al. 2001 p3|). With our definitions, those numbers, v and X, are twice the 
usual Ekman numbers, and they then quantify the ratios between the viscosities and the Coriolis force. 

- H and E respectively the dimensionless enthalpy and the dimensionless external accelerations, both scaled by 
the inverse of a. 

Concerning the latter quantities, we cut them in two parts and use for variables the difference between their present 
values and their values in the steady state. For the background parts, we then have 

So = - W + yli^— with VH Q = E . (7) 
2 a 

This equality enables the background to be solution of the NSE and can be solved separately. From now until the 
end, we consider that this condition is realized, and we forget about So and Hq (except in the mass conservation 
equation where the latter still appears as an external parameter) . 

Finally, some words about the Newtonian gravitational field. Inertial modes are current oscillations associated with 
small density variations. In this context, it is natural to use the Cowling approximation (Cowling 1941 p4|) which 
consists in forgetting fluctuations of the gravitational field (see the relevance of this approximation for Newtonian 
r-modes in Saio 1982 p5f). We will apply it for the non-linear study. Nevertheless, for the Newtonian linear study. 



we shall see later that depending on the way mass conservation is treated, it may be not necessary (see Section [II 
for more details). In this case, we have 

S = E*o + o with a = -X/5U (8) 
where SU is the Eulerian perturbation of the dimensionless gravitational potential. 

B. Dimensionless equations of motion for the spherical components of the velocity 

The equations of motion for the spherical components of the velocity in the orthonormal basis associated with 
(f , i?, ip) are given by 



d T V r + 8 V V r + [V • Vj V r - § (V$ + V%) - 2 sm r &V tp + d^h = 
% ((^a) pdtV r ) + s (AV r - $ (a^fc (smW) + <fad v V v + V r ) 



d T Vg + d^Vo + (V ■ V) Vg + | (V#V r - V* COttf) - 2 co* \ - [ -i)„h = 



2 ((d e a) (Vtf + | (d$V r - Vi)) + a (aVo + £ (28»V r - ^ - 3$$d v V v )) (9) 

+|^((¥ + l)v^))+^, 
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3 T V V + d v V v + (V-V)V (P + I (V v V r + V V V<> cot[0]) + 2 (sin <d V r + cosi?^) + ^^d^h 



where 



1 12 (t f\ i 1 o /„• .Ofl.M i ' ;)-' 



£ q £ 2 smtf f sm $ 

is the scalar Laplacian, while 

\ / t t sin 17 



and 



V • ?=ifl 4 {e V r ) + -r-?-^ (sini?^) + -J—d^. 
^ £ sin if 5 sin 17 



In the above equations, we assume s to only depend on £ in the (£, tf,^) system of coordinates. This assumption 
is linked with the slow rotation limit we should apply from now. This approximation supposes that the star is 
slowly rotating compared with its Kepler frequency. As explained in the introduction, observational data make this 
assumption credible, while known pulsars seem to rotate with a velocity smaller than a third of their Kepler velocity. 
In the lowest framework, only terms linear in fl are kept, and the background star is supposed to have preserved its 
spherical shape. Yet, to improve the slow rotation approximation, one can go a step further and take into account 
the deformation of the star, assuming it is a small order effect (or small quantum number / effect). To make this 
improvement, it is sufficient to introduce a new variable x[d that is equal to the unity at the surface, that coincides 
with isosurfaces of pressure, enthalpy and density, and that decomposes on m — Legendre functions 4 . With our 
algorithm, the easiest way to proceed would be to keep the spherical basis for the vectors, and to express all spatial 
operators depending on (£, <p) as the sum of the new operators depending on (y, ip) and of terms to interpret as 
applied accelerations. More details can be found in Bonazzola et al. 1997, 98, 99 |3(|, |37j and |Q. This will not be 
done in this article, for reasons that were also explained in the introduction. Finally, note that we do not really solve 
these equations. Indeed, we use the Helmholtz theorem (see for instance Morse and Feshbach 1953 (39)) that says 
that any vector of K 3 can be in an unique way written as the sum of a divergence free vector and of the gradient of a 
potential, given its normal component over the boundary: 

V = V-0 + V A V. (10) 

Instead of working with the above equations [Eq.(^)], we use the equation on the scalar potential tp and the equations 
on the 1/ and ip components of the divergence free vector V A V. Most of the time, these equations cannot be reached 
analytically. The way to proceed is then to separate the potential and the divergence free parts of the initial equations 
by numerically solving Poisson like equations obtained by taking their divergence. The reader can find in the appendix 
more details about these algorithms. 



C. Characteristic numbers 



To better understand what are the dominant processes in the dynamics of NS, it is worth estimating characteristic 
time scales and associated numbers. The acoustic time, i.e. the duration of the acoustic waves travel across the NS 
is by far the shortest. For a typical NS, it is about 2.10 -4 seconds 5 . As we are dealing with slowly rotating NS, 
the period of rotation should be more than ~ 2.10 -3 seconds (= 500Hz). It follows that this time is at least one 
order of magnitude greater than the acoustic time. It is then the same for the typical period of inertial modes, their 



4 Restricting to the case m = is sufficient for an isolated NS, but in a binary system, m = 2 should also be included. 

5 This time is of the same order of magnitude as the inverse of the frequency of the fundamental pressure mode, the so-called p-modc. 
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frequency being (at linear order) proportional to fl. 

Another time scale is the viscous damping time associated with viscosity: 

n a R 2 



T 



max [77, C] 



T 



2 max [E s ,E b ] 0, 



where E s and Eb are respectively the Ekman number for shear and bulk viscosity. In other words, the Ekman 
number must be interpreted as characteristic of the ratio between the period and the viscous time. This number is 
typically less than 10~ 7 and the viscous time is more than 7 orders of magnitude greater than the period (see Cutler 
et al. 1987 |4(j] for more precise values). A flow will be said "rotation dominated" if the above Ekman and Rossby 
numbers are small compared to the unity. Note that another usual hydrodynamical number appears with them: the 
Reynolds number, prodrome of turbulence. In a rotating fluid, it is defined by the ratio between Rossby and Ekman 
numbers, or in any fluid by the ratio between the non-linear and the viscous terms. 

The last typical time to evaluate here is the instability rising time T g associated with the RR force. To get an 
idea of its value, we come back to the dimensionless RR acceleration, formulae (0) and definitions (H). Analytical 
calculation with V 1 being the I — m linear r-mode (with time dependent amplitude) 

V = A[t}-^ k m f A VF m J (11) 



rn 



or in the spherical orthonormal basis 

V = A[t] 





^(sintf)" 1 " 1 sm[m(p] (12) 
£ m (sin$) m_1 cos[i9] cos[m<y9] 



(where K means the real part of the complex function), gives in dimensioned variables at the lowest order for the 
m = 2 r- modes 

n tr . GR 7 n 6 n 2 16 tt A , _ 1 A , , 

d tA [t] = g7 38527 ^]=^^]- ( 13 ) 

For a typical spherical NS with R = 10 km, M = 1.4 M Q and fi - (2tt)200 Hz, T g is something like 10 8 periods 
of the NS. Thus, depending on the viscosity given by the EOS, it will or will not be larger than T v , and then, the 
inertial instability will or will not be relevant. At this step, a new typical number seems natural to introduce. We 
shall propose to call "Chandra number" 6 the ratio between the viscous time T v and the rising time T g . In the same 
spirit as the Rossby and Ekman numbers, it should also quantify the ratio between the viscous and RR forces: 

3 8 5 2 7 c 7 max[77, (] ' 

Note finally that with this factor ahead of the physical parameters, we ensure the bifurcation value of Ch to be of 
the order of the unity, at least for the linear I — m — 2 r-mode. Indeed, this point is easily illustrated by looking at 
NSE for the linear I — m — 2 r-mode with time dependent amplitude and a shear viscosity of the form s — 1 — £ 2 . 
This viscosity that vanishes at the surface implies no need to add more boundary conditions (BC) and gives the exact 
(at the lowest order for the RR force) differential equation for the amplitude: 

9 t Alt] = ±-A [t] (l-A 



Chandrasekhar was the first who studied the gravitational radiation driven instability for the I = m = 2 fundamental modes of uniform 
density MacLaurin spheroid in 1970. See In]. 
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TABLE I: Characteristic numbers implied in the dynamics of inertial modes of rotating NS. Note that the definition of 
Chandrasekhar number is adapted to be of the order of 1 for linear I = m = 2 r-modes. 



Numbers 


Definition 


Analytical expression 


Ekman 


Ratio between P and T v or between the viscous term and the Cori- 
olis force 


E ~ v/ (nQR 2 ) 


Rossby 


Ratio between the typical velocities of the mode and of the back- 
ground fluid or between the non-linear term and the Coriolis force 


Ro ~ W/ (RQ.) 


Reynolds 


Ratio between the Rossby and Ekman numbers or between the non- 
linear and viscous terms 


Re ~ nWR/u 


Chandrasekhar 


Ratio between T v and T g or between the RR force and the viscous 
term 


Ch ~ n 2 GR 9 Q 6 / (c 7 u) 



TABLE II: Typical values of the different time scales implied in the dynamics of inertial modes of rotating NS. 



Time scale 


Definition 


Analytical expression 


Typical value (or range) 


Acoustic 


Travel of acoustic waves across the NS 


rp R 
a ~ 


~ 2.1(T 4 s 


Inertial 


NS's period (same order as the period of 


p _ 2vr 

r ~ n 


> 2.10" 3 s 




the linear inertial mode) 




Viscous 


Damping of inertial modes due to viscosity 


T v ~ 1/ (2EQ) 


> 10 7 P 


Gravitational 


Growing due to RR force 


T g ~ c 7 P/ (GMR 4 Q. 5 ) 


~ 10 8 P 



where one easily sees that with this shape, the viscosity wins the battle against RR force for Ch < 2. All the previous 
numbers are gathered in Table Q. 

To end with this short discussion, we should insist on the most important conclusion that is summarized in tab. 
the above figures show how stiff is the numerical problem of finding a dynamical solution of the NSE in this 
framework. This is a physical situation in which several very different time scales appear. In order to get an efficient 
and accurate code, some approximations have to be made. 



III. MASS CONSERVATION, BOUNDARY CONDITIONS AND APPROXIMATIONS 



To be consistent with the variables we chose in NSE, we have to write the mass conservation equation with the 
dimensionless enthalpy. With the decomposition explained at the end of the Section [I A , we have 



H 



(15) 



The first term Hq is the background, the second h corresponds to the mode itself and a still quantifies the non- 
linearity. This gives 



(d T +d v )h+(y- VH + T[n]H V • V ) + a (v ■ Vh + T[n]h V • V 







(16) 



where T[n] = 



ding 
d In n 



In the polytropic case, P — kri 1 , T is constant and reduces to T = 7 — 1. Exception done of the 



case 7 = 1 where the enthalpy is a logarithm. In the linear case, we have to neglect the last part of this equation, i.e. 
to do a = in Eq.((l|). We shall now describe some different choices that can be made to deal with this equation. 
The reader can find in appendix ^ the algorithms to implement all the following schemes in the framework of spectral 
methods. 



A. Solving the exact system of equation 



If one tried to solve numerically the exact Navier-Stokes and mass conservation equations, one would have two 
main possibilities. The first would be to employ an explicit scheme. But by this way, the Courant conditions for the 
following of the acoustic waves would impose a time step very small compared with the period of the star. This would 
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almost forbid any hope to make evolutions during durations long from the inertial modes point of view. The second 
way to proceed would be to use an implicit scheme for the divergence free part (see the appendix) . Indeed, this would 
make it possible to take a time step not too small compared with the period of the NS. But here the problem would 
be to estimate the errors done. Hence, in a first study, some approximations can be introduced to solve the problem 
in a easier way. 



B. The divergence free approximation 

To avoid the problem of solving acoustic waves for better studying inertial modes, the easiest solution is to use the 
divergence free approximation. It consists in replacing the usual mass conservation equation by V • V — 0. By this 
way the time step can be chosen larger than in the solving of the exact system and then allows the following of the 
inertial modes themselves. From a numerical point of view, this is a fast and robust approximation quite useful in the 
exploring phase of the numerical work. Yet, it should be noticed that this drastic approximation may be not too bad 
for inertial modes in NS. Indeed, the I — in = 2 r-modes which are the most interesting for GW have a divergence 
free limit in the linear order. Furthermore, the latest figures about bulk viscosity (Haensel et al. 2001 153|) that are 
quite huge, could mean a fast damping of all modes, except of those that are divergence free. Finally, note that in 
the linear limit of EE, this approximation gives an evolution of the mode V that is independent of the background 
star if it is rigidly rotating. 



C. The anelastic approximation 

Yet, instead of being quite "savage" with the equation and imposing the divergence free condition, one could look 
for a cleverer way of doing physics and think about inertial modes. A main feature of these modes is that their 
frequency is of the same order as the angular ve locit y of the star in which they occur. As their damping and growing 



times are larger than their period (c/. Section II C), one would like to take it for a characteristic time scale. From 
practical and numerical points of view, it means a time unit (or time step) choice not too small compared with the 
period, in order not to waste memory and computational time calculating non relevant physics. From a physical angle, 
the philosophy is to neglect acoustic waves by assuming that time derivatives of the pressure and density pertur- 
bations do not play a key role in the phenomenon. It can be done in a consistent way with the anelastic approximation. 

The anelastic approximation was first introduced in atmos phe ric physics by Batchelor in 1953 p| and then derived 
from a rigorous scale analysis by Ogura and Phillips in 1962 |42| who gave it its name. In astrophysics, it appeared in 
1976 in a paper by Latour et al. p3[ concerning convection and was then widely used in that field and others (such 
as stellar oscillations) where one can neglect temporal variations of the perturbation in density but not necessarily 
spatial ones. For a recent critical approach of this approximation within astrophysics, one can read Dintrans and 
Rieutord 2001 @ and Rieutord and Dintrans 2002 @|. 

In the Eq.([l^), anelastic approximation consists in neglecting (d T + d v ) h which is the time derivative of the enthalpy 
in the rotating frame. By this way, we cut acoustic waves and then have 

V ■ VH + T[n]H V • V + a (v-fh + T[n]h V ■ V ) = 0. (17) 

As a conclusion on this approximation, we would like to insist on a particularity of the linear case. With anelastic 
approximation and linearization of all equations, we obtain an equation for the mass conservation that does not 
depend on the Eulerian perturbation of the enthalpy But, taking the curl of the linear EE (or NSE) equation gives an 
equation with exactly the same feature (remember that the curl of a gradient is zero) with an interesting additional 
fact: it neither depends on the Eulerian perturbation of the gravitational potential. Furthermore, it is easy to verify 
that for a background enthalpy depending only on £, the boundary condition is that the radial velocity should vanish 
at the surface. Thus, the situation is that finding the Eulerian perturbation of velocity should give exactly the same 
result, whatever the hypothesis on the Eulerian perturbation of the gravitational potential. It means, we do not need 
Cowling approximation with anelastic approximation in the Newtonian linear case. Cowling approximation would 
play a role only if we wanted to find what is the enthalpy and what is the gravitational potential in the source term 
of the gradient part of EE. 
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D. Boundary conditions 



The only missing information is now the boundary conditions we chose. In an actual NS, the fluid core is supposed 
to be surrounded by a more or less rigid crust, an ocean and an atmosphere. Their physics is by itself quite a complex 
subject (see for example Haensel 2001b But even for a simple toy-model, for instance a crust made of only 

one type of nuclei, depending on the temperature and on the amplitude of the motion of the inner fluid, the physical 
state of the crust can be quite complex to describe, something like icepack on the sea (see for instance Lindblom et 
al. 2000 0]). There is no need to explain how difficult would be to translate this BC in a mathematical language... 
This is the reason why, for the EE, we then chose to begin with two different and extreme BC to get an idea of the 
limit cases. The first is the free surface BC, i.e. the absence of any crust. The second is the presence of a rigid crust 
at £ = 1, if we take the inner radius of the crust for the typical length R. In the first case, one has H |g=i = and 
in the second V r |f=i = 0. The numerical w ays to take into account those BC can be found in the appendix B3. For 
the NSE, as it was already said in Section HC, in order to avoid the need of more BC, we chose a degenerate shear 
viscosity of the form s = 1 — £ 2 . 



IV. TEST AND CALIBRATION OF THE CODE 



Now that we have presented our physical framework, we will focus on the code. As it was already mentioned in 
the previous sections, it uses spherical coordinates and spectral methods to solve NSE. More precisely, it solves the 
equations coming from the decomposition of NSE into a potential and a divergence free parts. We will not give 
more details here and send the reader to the appendix for informations about the algorithms and spectral methods. 
The rest of this article, and the discussion about numerical stability in the appendix, are devoted to the linear study. 
Non-linear work is still in progress and will be described later. Furthermore, in this section, we only deal with rigid 
rotation in order to have some analytical solutions to make tests with. 



A. Conservation of the energy 

The first test of the code was obviously the free evolution of the linear I = m r-mode in a rigidly rotating inviscid 
and incompressible fluid. The vector defined in Eq.(|l2]) is indeed an eigenvector of EE whose frequency in the inertial 
frame is = — ^ m ~^™ +2 ^ ■ Moreover, note that the absence of radial component implies that both free surface and 
spherical rigid crust BC are automatically satisfied. Yet, for all the preliminary calculations using the divergence free 
approximation, the code was built to work with the rigid crust BC. 

The great advantage of spectral methods is to change any linear spatial operation in linear algebra calculations, 
which can be done exactly. As the velocity of Eq. (|lj) has a very small number of coefficients in the reciprocal space, 
it is easy to verify (for example by directly looking at the time evolution of those coefficients) that there are only 
errors due to round-off and to time discretization. In figure ([!]) is illustrated the time evolution of V& at the equator 
for the m = 2 linear r-mode. The spatial lattice was of the shape (8x6x4) for (N r , N#, N v ) with 100 time steps 
per period of the mode. The duration of this run was chosen for the evolution to last exactly 100 times the expected 
period. This calculation was done on a DEC Alpha Station with 500 MHz processor and took 217 seconds of CPU 
time. Comparison of the amplitude at the first step and at the last steps shows the growth of the amplitude due to 
errors that come from the time discretization. This is easier to see in figure (Q) where the time evolution of the error 
in energy appears. The different curves illustrate results obtained with the same grid and duration, but in runs with 
different numbers of time steps per oscillation. Here are pictured results for 100, 150 and 200 [see also associated 
power spectra in figure (g)]. The last two calculations on the same computer took respectively 269 and 342 seconds 
of CPU time. 

In figure (Q), we drew for several runs the relative error in the energy per oscillation of the mode, versus the 
number of time steps per oscillation. Both are in logarithmic scales. This error appears to exactly varies as At 3 and 
then as the inverse of the cube of the number of steps per period (regression slope of —3.01). This is due to the 
second order scheme we use to solve the NSE or EE. Indeed, when the energy is calculated, the error of order At 2 



In fact, symmetries of the spherical harmonics are used and the effective value of l max is roughly twice the number of points in i5. 
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FIG. 1: Time evolution of the i? component of velocity on the equator for the linear m = 2 r-mode in a incompressible 
and rigidly rotating fluid. The run is supposed to last exactly 100 times the period of this mode which is The beating 
phenomenon which can be seen is due to the resolution of the graphical tool. See the power spectrum in figure . 



done on the velocity is multiplied by a source term linear in At. 

As we are now only dealing with the linear code, we will not discuss the stability of the non-linear code, nor give 
test pictures corresponding to modes with azimuthal numbers different from 2. Indeed, there is a coupling between 
different values of m only in the non-linear versions of NSE and EE. The modes that are the most interesting for 
GW are m = 2 modes and in the following, we will always choose this kind of initial data. Nevertheless, the same 
calibration tests for the other linear I = m r-modes were done and gave the same power law relations between and 
the number of steps per oscillation. Finally, we also verified that the error in the phase of the modes is proportional 
to the square of the number of steps per period. The relative errors in the phase for the m = 2 mode with 100, 
150 and 200 steps per oscillation were respectively 1.64 10~ 3 , 7.33 10~ 4 and 4.12 10 -4 . It is easy to verify that the 
logarithms of these numbers are on a straight line with a slope very close of 2. 

We shall see now how we tried to improve the conservation of energy. For the free r-mode, the energy is expected 
to be exactly conserved. Furthermore, we know our numerical error in the velocity is of the order of At 2 . Remember 
that it comes from the second order scheme used to calculate source terms: 

gj+l/l = 35 J - SI' 1 (18) 

The basic idea is thus to modify this quantity, in the coefficients space, with some additional term of the order of 
At 2 , in such a way to retrieve the conservation of the energy without changing the error in the velocity. We then 
modified the Eq.(n8|) and took 



gj+1/2 = (3 + g )SJ-(l + e)5J- 1 



(19) 



where e is of the order of At 2 . It is calculated using values of the energy at the last two instants and in such a way to 
"impose" on the energy to be conserved. Obviously, when there are forces such as viscosity or RR, their power has to 
be taken into account in the energy balance equation. In fig.(^|) are the same curves as in fig.(||) but in logarithmic 
scales. We added the result of this "improved conservation of energy" for the run with 100 steps per oscillation. As 
the correction is local in time (done at almost each time step), it leads to a time independent error that should be 
able to remain the same as long as one could wish. This calculation took 320 seconds of CPU time. 
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5E/E versus t for 100 oscillations (without any correction) 
m-2 r-mode with grid 8x6x4 



! I 1 I 

100 steps per oscillation 

150 steps per oscillation 

200 steps per oscillation 




FIG. 2: Time evolution of the relative error in the energy before any improvement of the conservation of energy. The different 
straight lines correspond to runs with the same duration (100 periods of the m = 2 linear r-mode), same spatial grid (8x6x4 
points) but different time steps. We see the linear variation of these errors with time. It shows that exception done of round-off 
errors (which are so small that they do not appear here), the main error is deterministic and then well controlled. 



Power spectra of V for 100 oscillations (without any correction) 

m-2 r-mode wilh grid 8x6x4 



1 1 1 1 1 1 1 

100 steps per oscillation 


1 1 1 1 1 1 1 1 1 1 1 


150 steps per oscillation 




200 steps per oscillation 
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FIG. 3: Power spectra of component drawn in figure (|lj) and of the same component calculated with 150 or 200 steps per 
oscillation. See also figure (Eh. 



B. Modes driven to instability in the rigidly rotating case 



After trying to improve the conservation of energy for the free case, we did tests with a RR force included. For 
instance, we switched on the RR force in order to drive toward instability the linear m = 2 r-mode. As it has already 
been implied in Section II C, we did it by using a modified version of the formulae ^ (see this section). For these 
calculations, we compared numerical results with analytical calculations, both reached with the same approximation. 
We took for the background star a homogeneous NS with M = 1.4M Q , R = 10 km and £1 ~ (2tt) 500-ffz. In hg.(|^) 
appear two time evolutions of the ratio between the energy and the inertial energy for 100 oscillations of the mode. 
The first calculation was done without the enhanced conservation of energy and the second with this improvement. 
The analytical calculation gives for the ratio a final value of the order of 1.0575 whose difference with the improved 
numerical value, 1.0582, is less than 10~ 3 . 
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FIG. 4: On this picture clearly appears the fact the error exactly comes from the adopted numerical scheme with an error that 
vanishes as the third power of the time step. The relative error and the number of steps per oscillation are both in logarithmic 
scales and the regression slope is then —3.01 ± 0.06%. 




FIG. 5: Here are pictured the same curves as in fig.(|2j), but in logarithmic scales. We added the curve corresponding to a 
run with 100 steps per period of the m = 2 linear r-mode with the correction of the conservation of energy described toward 



Eq.@ 



The next step was to switch on the RR force, but without the linear r-mode for initial data, even if we still assumed 
that the fluid was divergence free. Instead, we chose for the initial velocity a random Gaussian noise with a first 
moment equal to and a second moment equal to 1. It should be noticed that in this case and others where the initial 
energy is not only distributed within quite low frequency modes, the above method to improve the conservation of 
energy is no longer appropriate and we shall not use it. Moreover, even if this method does not change the angular 
momentum when applied to the m = 2 part of the velocity, it can if one naively applied it to the m = part. But, 
as we already mentioned, here we only deal with the linear code and mainly with m = 2 velocities. And in this case, 
we exactly control the error done, without changing the angular momentum. We drew in figures (|) and © the 
time evolutions of the r and i) components of the velocity and their power spectra. In fig. (^|) is illustrated the time 
evolution of one of the two independent components of the S%j tensor, with the associated power spectrum. These 
calculations were done with a lattice of the shape (12 x 8 x 4), with 200 steps per oscillation and lasted 200 times 
the period of the m = 2 r-mode. The included RR force corresponds to what should exist in a NS with M = 1.4M©, 
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Time evolution of energy with reaction force 



Without correction 
With correction 




FIG. 6: Time evolution of the ratio between energy and initial energy in an inviscid and incompressible rigidly rotating fluid 
with a RR force corresponding to what should exist in a NS with M = 1.4M©, R — 10 km and Q ~ (2ir) 500Hz. The first 
curve corresponds to the basic code and the second one to the code with improved conservation of energy. The expected final 
value is about 1.057. The difference between this analytical calculation and the corrected numerical result is less than 10~ 3 . 
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FIG. 7: Time evolution of the radial component of velocity in an inviscid and incompressible rigidly rotating fluid with 
Gaussian noise for initial data and associated power spectrum. A huge RR force acts, but as expected, no polar counter part 
of the mode grows. 



R = 10 km and £1 ~ (2ir) 1050Hz. Thus, this calculation is not really physical (due to the use of the slow rotation 
approximation) but just aims at giving a qualitative idea of what should happen with physical values. As expected, 
there is an axial mode (no radial velocity) driven to instability with exactly the frequency (w — |fi) and coefficients 
(I = 2) of the r-mode. Moreover, a single look at the figure (|^) shows that even when the velocity is mainly noisy, the 
tensor that plays the key-role in the instability is quite smooth. In this figure, we also drew a zoom corresponding to 
tfi < 50 and the associated power spectrum. It shows two frequencies in this component of the tensor. One of them 
is the unstable r-mode with 3w/f2 = 4 and the other (with 3w/tt ~ 2.70) disappears with a longer run, as the scale 
is adapted to the growing mode. Yet, even in the spectrum of the full evolution, a trace of it can still be seen. We 
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FIG. 8: Time evolution of the i? component of velocity on the equator in an inviscid and incompressible rigidly rotating 
fluid and associated power spectrum with same initial data in fig- (0) - Here, the I — m = 2 linear r-mode is clearly driven to 
instability by the RR force. As expected, its frequency is w = Others peaks are just marks of the initial noise. 




FIG. 9: Time evolution of one of the two independent components of the Sij\t] tensor that appears in the RR force. This 
calculation was done during the same run as the results in fig. (0) an d (§)• We can also see the almost monochromatic associated 
spectrum. Nevertheless, there is another frequency in this tensor that is not driven to instability. It can be seen more easily in 
the zoom of the time evolution or in the corresponding power spectrum. 



achieved exactly the same features with others noisy initial conditions, for instance a Dirac kick into the NS 8 . We 
also did the same calculations with other spatial lattices (up to 64 x 48 x 4), and this result did not change at all. 



We mean an initial velocity equal to anywhere except in an arbitrary point. 
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C. Test of the anelastic approximation 



The last calculations we did in a rigid and Newtonian background were to test the effect of the anelastic 
approximation. The idea was to compare the previous results obtained with the divergence free approximation and 
a rigid crust BC with results coming from the same initial conditions but with the anelastic approximation and the 
free surface BC. For the rest of this section, the background star is a polytrope with 7 = 2. 



First, we took for initial data the linear I = m = 2 r-mode and let it freely evolve. Here again, the BC do not 
play any role. The spatial lattice was (8x6x4) with 100 time steps per period of the mode. The time evolution of 
V/} at the equator [same as in fig. ([!])] did not show any difference with the divergence free case, even for the power 
spectrum. But this is quite obvious to understand when looking at both the Eq.(|l7|) and the EE. Indeed, we see 
that the linear r-mode, which is divergence free and has no radial component, is still an eigenvector of this system 
of equations. Then taking it for initial data in the divergence free case or in the anelastic approximation should give 
exactly the same evolution. We verified [see fig. (10)] that the radial component of the velocity does not grow and 
neither does its divergence. This calculation took 344.6 of CPU time without any optimization and with a very basic 
solving of the anelastic equation. 



The next step was the driving toward instability of a mode with noise for initial data. We took exactly the same 
conditions (duration, lattice, initial data) as in the divergence free case, exception done of the boundary conditions. 
For these calculations and all that follows with the anelastic approximation, we will use the free surface condition that 
is automatically satisfied (see the appendix for more details). The figures ([[l]), (|l2|) and (jl^) that show respectively 
the time evolutions of V r , of Vg and of a component of Sy give a kind of summary of the results: 

- the radial velocity still remains noisy and does not grow; 

- the $ component grows in the same way as in the divergence free case. There is a difference between the final 
amplitudes that comes from the fact that the coupling to the RR force depends on an integral on the whole star 
of a function that is proportional to the density (and then depends on the EOS). Note that the power spectra 
have the same secondary peaks as in the divergence free case. 

- the spectrum of the tensor shows that there is only one unstable mode which has again the frequency of the 
linear r-modc. 

The conclusion is then that the anelastic approximation gives results very close of those obtained in the divergence 
free case. This is due to the fact that, in spite of the presence of all modes in the initial conditions, the only growing 
mode is the r-mode that has no radial velocity and is divergence free. 



V. DIFFERENTIAL ROTATION 



As we already mentioned it in the introduction, Kojima first noticed (see Kojima 1998 p5| ) that relativistic r- modes 
may be quite different from what they are in a Newtonian background due to frame dragging. Yet, since the main 
reason for this difference is the modification of the NSE and the possible appearance of what is called a continuous 



spectrum (see Ruoff et al. 2001 |26 , ETj] and Beyer et al. 1999 [|48|), the same may append even in the Newtonian 
case if the star is not rigidly but differentially rotating. And there are several reasons why a NS may not be in rigid 
rotation. First, the birth conditions of the NS themselves. Secondly, the non-linear coupling of modes. Then, a 
possible drift induced by the existence of a magnetic field. As we are here only dealing with linear hydrodynamics and 
tests of the code, we will not giv e more details about those processes (and send the reader to the follo wing articles 
for more details: Spruit 1999 @, Rezzolla et al. 2000 @, 2001 [fl] and ||, and Schenk et al. 2002 §f|). What 
we have done is just to take as given that the background star is differentially rotating with quite an arbitrary law 
and to look at the influence of this law on the existence of the modes. Nevertheless, instead of asking the question 
"Is there any r-mode left in a differentially rotating NS?" that is not well defined, we decided to try to answer to two 
different and more precise questions: 

- is there anything growing when a RR force is applied on noise in a differentially rotating background? 

- what does happen to the linear r-modes if they are chosen to be the initial data in such kind of background? 



Some lights on these questions are in the following subsections, but we shall begin with some words about the 
modifications implied on the equations by differential rotation. 
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FIG. 10: Time evolution of the radial component of the velocity in the equatorial plane in £ = 0.5 with the linear m — 2 
r-mode for initial data. The background star is a rigidly rotating 7 = 2 polytrope and the calculation is done with the anelastic 
approximation and without any RR force. The run lasts 100 times the period of the r-mode. There should be no radial velocity 
(c/. text) and we see it comes only from round-off errors (double precision calculation with initial values of the order of the 
unity) . 




FIG. 11: Time evolution of the radial component of velocity with Gaussian noise for initial data and its power spectrum. The 
background is a polytrope with 7 = 2 and the anelastic approximation is used. As in fig.(^) a huge RR force acts, but as the 
star is rigidly rotating no polar counter part of the mode grows. 



A. Modifications 



Assuming differential rotation of the background star involves slight modifications of what we said in the Section |J. 
First, terms coming from the spatial derivatives of the non constant f2[£, 1?, <p] must be added in the NSE. Secondly, 
we shall now specify what exactly means f2 in the definition of dimensionless variables [cf. Eq.(|^)]. 

Concerning f2, we chose to take for a time unit the inverse of its value at the equator. This is uniquely determinate 
due to the fact we have always assumed that the rotation law is of the form Q[r, <p] = fi[r sin[i?]] or of the form 
f2[r, $, (p\ — Q[r]. The first case corresponds to what must be this law for the background to be stable with respect to 
the Newtonian EE, and the second case is in a way inspired by GR even if it is not a solution of the full Newtonian 
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FIG. 12: Time evolution of the i? component of velocity and associated power spectrum for the same calculation as in fig. 
([[!]). The only difference with the divergence free case [cf. fig.([s|)] is the final amplitude. Yet, this is not due to the anelastic 
approximation but to the effect of the EOS on the coupling coefficient of the RR force. 




FIG. 13: Time evolution of one of the two independent components of the Sij[t] corresponding to the same run as the results 
in fig.(|ll]) and (|l2|). Here can really be seen the fact that there is only one unstable mode with the frequency of the linear 
r-mode and that the anelastic approximation does not change this feature. 



EE. For more details see Section VI 



In the dimcnsionless NSE, the modifications induced by differential rotation are quite simple. First, d v is replaced 
with Q[r, where Q[r, i9] is the dimcnsionless profile of rotation. Then, we have to add new terms coming from 

e v r sin[#] V ■ V (&[r, i?]J that are in the spherical orthonormal basis 






V r rsin[&]d r n + V» sin[&]d#£l. 



(20) 
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B. Noise with huge RR force 



Once again, our goal was to stay as close as possible of the basic and well understood situation to minimize the 
number of unknowns. Thus, we took noisy initial data, put it in a differentially rotating background, switched on the 
RR force and looked to what was to happen. By this way, the question was not to look for the existence of r-modes, 
but simply to look for the existence of modes driven to instability by the RR force in a non rigidly rotating background. 

The only difference with the basic study done in the case of rigid rotation was that we had to choose a law for O. 
The first choice was very simple and corresponded to a background stable with respect to the Newtonian EE. It was 
of the form 

S7M] = W (l + /?„ r 2 sin[tf] 2 ) (21) 

where (3 n was a constant depending on the run and W calculated to have f2 [l, ^] =1. As explained above, we also 
tried a law inspired by GR: 

n[r,0] = W(l + (3 gr r 2 ) (22) 
or laws coming from the already quoted article by Karino et al. Q: 

"Ml= i ■ 7j£h- n 2 (23) 
H sink? J + Pl 



or 



ftM = - W f v xl/2 . (24) 



(r 2 sin[z9] 2 + /3, 2 



Since in all these laws, none of the free parameters [3 corresponds to a physical variable, we will not give 
quantitative results. It will be done in the relativistic study (and then in another article) where there is a parameter 
that physically quantifies the way the equations are far from the Newtonian EE: the compactness of the star. Here 
we will only discuss in a qualitative way the results obtained in all the previous cases and give a very representative 
example: what happens in a 7 = 2 polytrope with anelastic approximation (free surface) and the rotation law given 
by the Eq. (pl|) with (3 n — 0.4. Exception done of these choices, the calculation was done with exactly the same 
conditions as in the last study with noisy initial data in the previous section. 



The main difference with the case of rigid rotation is that the mode that grows is no longer purely axial. Indeed, 
we can see on fig.(|l4|) that the radial velocity is now also driven to instability. Moreover, comparing power spectra in 
this figure and in fig.([l5|) and ( |l6| ) shows that this is really a single mode with frequency very close of the frequency 
of the r-mode. By a single mode, we mean that this is exactly the same frequency for several physical quantities and 
several positions in the star. This is not a sufficient condition to talk about "discrete spectrum" , but it is a necessary 
condition. This unavoidable existence of a polar part of the velocity in a differentially rotating Newtonian star with 
barotropic EOS should be compared with the results achieved by Lockitch et al. Q in the relativistic framework. 
Indeed, in GR, the main reason for the coupling between axial and polar parts of the velocity is the frame dragging 
that is imitated in a Newtonian framework by differential rotation. Finally, note that the value of the frequency in 
the Newtonian case is depending on the way we chose to normalize fl. Then, to summarize all our calculation, we 
can say that we observed that the polar counter part of the mode appears as soon as there is differential rotation, and 
whatever the chosen law. Yet, the more the law for £1 is far from the rigid case, or in more pragmatic way the greater 
the free parameter is, the more the unstable mode has a polar counter part. 



C. Free evolution 



The other question we asked was: "What does happen to linear r-modes if they are put into a differentially rotating 
background?" . The idea in taking such kind of initial data was to have something quite close of an eigenvector, 
assuming if there is one it should be quite similar to the linear r-modc. 

Once again, we did not make quantitative calculations and postponed it to the GR case. As in the case of noisy 
initial data driven toward instability by a RR force, the free evolution of the linear r-mode showed the growth of 
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FIG. 14: Time evolution ol the radial component of velocity in a 7 = 2 polytrope with anelastic approximation and Gaussian 
noise for initial data. A huge RR force acts (something like what should exist in a star with angular velocity equal to 1000 Hz) 
and the background star is assumed to be differentially rotating with the law corresponding to f3 n = 0.4. We see that a polar 
counter part of the mode now grows. 




FIG. 15: Time evolution of the ■& component of velocity in the same calculation as in fig.(y_4J). The associated power spectrum 
shows that there is one mode driven to instability that is the same as in fig.(^4|) and has a frequency very close of the expected 
frequency of the r-mode. See also the fig.dl^). 



a polar part of the velocity. It can be seen on figures ([i~7|) and ( |l8| ) that respectively show the time evolution of 
the ratio between the polar energy and the total energy and the time evolution of the radial velocity at the point 
of coordinates (£ = i, , *9 = : |,(^ = 0). This results comes from a calculation done with a spatial lattice of the shape 
(64,48,4), the anelastic approximation (with a free surface) and a rotation law given by th e Eq.(|2l|) with (3 n — 0.4. 
For reasons that will be explained later, we included degenerate viscosity (c/. Section HID) with an Ekman number 
E a = 5. 10~ 5 in order to regularize the solution. The evolution was done to last 30 periods of the linear r-mode with 
100 steps per oscillation. In the figure (0), we see that, after a while, the ratio between the energy in the polar part 
of the mode and the total energy reaches a kind of stationary state with a coupling between different modes. The 
existence of this "hybrid final state" was verified during other runs with other physical conditions and is once again 
to compare with results achieved in GR by Lockitch et al. . 
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FIG. 16: Time evolution of one of the two independent components of the Sij[t] tensor that appears in the RR force. This 
calculation was done during the same run as the results in fig.(|l4|) and (^). We can see the almost monochromatic associated 
spectrum with the same frequency as in the previous spectra of the unstable mode. 



Concerning the existence of modes, looking simultaneously at fig. (|l8|) , (19) and ( p0[ ) shows that apparently one 
single mode mainly appears in the reaction force (or in the component of the corresponding tensor) even if both axial 
and polar part of the velocity contain several modes. Yet, the spectrum of the Sij tensor is quite noisy due to the 
fact that the RR force is not switched on and that the run is short. We verified this feature during other calculations. 
The conclusion is that for small values of the (3 parameter (these values are depending of the chosen rotation law), 
the main effect of differential rotation on a "free linear r-mode" is to give it a polar counter part and to widen its 
spectrum. But for larger values of the parameter, even if the velocity's spectrum becomes very noisy, the Sy tensor 
is always less noisy. Yet, here we will not give more details about this points, the linear Newtonian study being not 
really interesting from the NS point of view. 



To end with this section, we will mention a quite amazing feature found in all our calculations and illustrate it 
with an example coming from the previous free run of a linear r-mode put in the differentially rotating background. 
We always found that quite fast (in something like 15 periods of the linear r-mode) the main part of the velocity 
is most of the time concentrated in a region close of the surface of the star and of the equator. This is illustrated 
in the figure ( pl| ) where we drew the shape of the (f component of the velocity versus the radius of the star for 
several values of the $ angle. Note that we got the same results with different BC and even when we add viscosity 
to regularize the solution (and this is the case in this calculation) and to be sure this is not a numerical artifact. In 
fig. ([22]) appears the shape of $ component for an evolution done with the same grid, viscosity and initial data but 
with f3 n = 0.8. With such a huge value of the parameter that governs the rotation law, there is a lot of different 
modes that appear in the velocity spectrum. But what seems to be the most important result is the very high 
concentration of the velocity near the surface. This is analogous to the results achieved by Karino et al. in 2001 [[[J. 
In the conclusion, we will shortly discuss what is, in our point of view, a possible important repercussion of this result. 



VI. STRONG COWLING APPROXIMATION IN GR 



NS are among the most relativistic macroscopic objects made of usual matter in the Universe. Moreover, 
the instability of r-mode is due to GR. It is then quite natural to try to study r-modes in the GR frame- 
work. The problem is that it is quite difficult to evolve dynamically the full GR and hydrodynamics equations 
taking into account GW. This is the reason why we decided to use the strong Cowling approximation at the beginning. 

Here we will only give some preliminary results obtained in the GR case. A more complete and physical study 
will be in a following article in preparation. Our goal is just to show that the above conclusions concerning the 
appearance of a polar part of the velocity and the "concentration of motion" for Newtonian differentially rotating 
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Fraction of polar energy versus time for 30 oscillations 
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FIG. 17: Time evolution ol the ratio between energy contained in the polar part of the velocity and the total energy of the 
mode in the anelastic approximation with a free surface. The background star is assumed to be a 7 = 2 polytrope differentially 
rotating with the law corresponding to (3 n = 0.4. A viscous term coming from the degenerate viscosity with E a — 5. 10 -5 
is included. At the beginning of the evolution, there is energy coming from the purely axial initial data (the linear r-mode) 
and going to the polar part. Then this energy oscillates from the polar part of the velocity to the axial part and back with a 
constant mean value. 




FIG. 18: Time evolution of the radial component of velocity with the anelastic approximation and a free surface. This 
calculation was done during the same run as in ng.(|l7|). The initial data are the linear r-mode that is no longer an eigenvector 
and radial velocity quickly appears. Several different modes are present. 



NS still hold in the framework of relativistic rigidly rotating NS. 

As we are studying slowly rotating NS, which remain spherical, a very convenient and good approximation is to 
assume that the 3-space is conformally flat (see Cook et al. 1996 p4fl). Furthermore, as we are using the strong 
Cowling approximation, there is no GW even without this approximation. Then, in unit c = 1, the infinitesimal 
space-time interval can be written 

ds 2 = - (N[r} 2 -a[r] 2 r 2 sink?] 2 A^[r] 2 ) dt 2 - 2a[r} 2 r 2 sink?] 2 N*[r]dt dip + a[rf dr 2 

where, following the 3 + 1 formulation, N[r] is the lapse function and N v [r] the third contravariant component on the 
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FIG. 19: Time evolution of the i? component of velocity in the same calculation as in fig.(|18[). The curve and the associated 
power spectrum show that there are several modes, but only one of them appears in the Sij tensor with quite a huge amplitude. 
Moreover, it has a frequency very close of the frequency of the linear r-mode. See also the fig.(^o|). Note finally that the same 
frequency also appears in the spectrum of the radial velocity. 




FIG. 20: Time evolution of one of the two independent components of the Sij[t] tensor that appears in the RR force. This 
calculation was done during the same run as the results in fig.(|l8|) and (jis|) . We can see there seems to be one main frequency 
and a second one of smaller importance. The first one has exactly the same frequency as the unstable mode that appears in 
the previous figures and calculations where the RR force was on. 



spherical coordinate basis of the shift vector, a[r] is the conformal factor and finally df 2 is the infinitesimal interval 
in the 3-space. In this equation, we use the isotropic gauge and the slow rotation limit to assume that all these 
unknown functions are only depending on r (and this is the reason why, in the Newtonian study, we called "inspired 
by GR" the law of differential rotation in which the functions are only depending on r). In Hartle 1967 |5jJ, the 
system of coordinates is not the same as the one we use, but it is easy to show that the conclusions are identical since 
the mapping between the two systems is quite trivial. Furthermore, from the same article we know that TV^r] is of 
the first order in f2. 



To get the equations of motion, the relativistic equivalent of EE, we write the energy tensor of a perfect fluid 

= {p + p)Um v + pg>" / (25) 
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FIG. 21: <p component of the velocity versus the radius for several values of We can see the main part of the velocity is 
concentrate in a region near the surface and close of the equator. This "snapshot" of the profile was taken at a moment when 
the derivative versus the radial coordinate of the velocity is quite huge at the surface. 
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FIG. 22: •& component of the velocity versus the radius for several values of Here the concentration of the motion is higher 
than in the previous calculation. (3 n was 0.4 in fig.(^l|) and is now 0.8. Moreover, this calculation lasted three time longer. 
Nevertheless, note that here we did not choose to draw the profile when the derivative is huge and just drew it at the final 
instant. 



where p is the energy density, p the pressure and f/* 4 the 4- velocity of the perfect fluid. Then, we apply the conservation 
of energy: 

V^T*" = 0. (26) 

If we define the generalized enthalpy: 

dP 

dH= — , (27) 
P + P 

it gives 

U^XTVyH + Vff + WVvU 11 = 0. (28) 
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It is well-known that due to thermodynamics, the system of equations obtained by adding the baryonic number 
conservation V M (n Z7 M ) (where n is the baryonic density) to those four equations is a degenerate system. Then, we 



chose to work with the baryonic number conservation and the projections on the 3-space of Eq.(|28| 



Following the Newtonian case, we just have to add Eulerian perturbations to the rigid rotation and to linearize 
the equations with respect to both the amplitude of these perturbation and f2. For the generalized enthalpy, the 
definition of the perturbation is obvious, but for the full 4- velocity, we use the well-known results about rigid rotation 
of relativistic stars and write it as 



N[r 



1 



SU°[t,r,i9,(j)} 
5U r [t,r,$,, 

N[r] \ 2 SU^ [t,r, ■&,<!>] 



(29) 



Note that in this equation, 5U° is not a dynamical variable but is determined according to the constraint that 
is a 4-velocity (U^U^ = —1). Furthermore, SU r ,SU and 5U V are not contravariant components of a 4-vector but 
convenient variables that are the components of a 3-vector on the orthonormal basis associated with the spherical 
system of coordinates for the flat 3-space. It enables us to write the motion equations in a way very similar to the 
Newtonian EE. Indeed, writing the 3-velocity on the orthonormal basis associated with the spherical coordinates as 



V= 



su-° 



(30) 



and defining on the same basis 



(Q - A^[r])cos[tf] 
#[t,r,i?, <f>}= { - sin[i9] (O - Nf[r] + A[r]) 




(31) 



with A[r}= (n - N^[r}) 



d hi 



dln[r] 



rN* [r] 
2 



where ' is the derivative versus the radial coordinate, we have 



(dt + ndp) V+ 2zuAV + Vh = 0. 



(32) 



This equation is very close of the Newtonian EE, the main difference being that the 3-vector that appears instead of 
n is now depending on the coordinates (as in the case of differential rotation) but no longer parallel to the rotation axis. 



Concerning the baryonic number conservation, writing it in a Newtonian like way, we have in the slow rotation 
limit 

(d t + nd v ) n+ f^jh div(nV) =0 (33) 

where n is defined as n=n A^ 2 a, n being the baryonic number density. The natural generalization of the anelastic 
approximation is then 

div fiv) = 0. (34) 

As we are using the strong Cowling approximation, the background star appears in the equations of motion only as 
"external" (from the point of view of the mode) data. For any relativistic calculation, what we do is to calculate the 
background configuration using the already existing code illustrated in Bonazzola et al. 1993 Q and then to use the 
resulting lapse, shift, conformal factor and their derivatives with respect to the radial coordinate in our equations. 



Here we will only give one single example, showing we have in this relativistic case with the strong Cowling 
approximation results similar to those obtained in the case of the Newtonian differential rotation. We took a 7 = 2 
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Fraction of polar energy versus time for 30 oscillations 




FIG. 23: Time evolution of the ratio between energy contained in the polar part of the velocity and the total energy of the 
mode with the strong Cowling and anelastic approximations and the free surface BC. The background star is a 7 = 2 relativistic 
rigidly rotating polytrope with 1.74 solar masses and a radius equal to 12.37 km. The ratio between its kinetic energy and its 
mass energy was HT 16 . We see that in spite the fact the initial data are the linear Newtonian r-mode that is purely axial, as 
soon as the calculation begins, a polar counter part of the mode appears as implied by Lockitch et al. |2|. 




FIG. 24: Time evolution of the radial component of velocity at the point (5, -jjO) for the same calculation as in fig.([23|). The 
initial data are the linear r-mode that is no longer an eigenvector and radial velocity appears very fast. The power spectrum 
of the time evolution shows several peaks, but one of them has the same frequency as the modes of fig.(p5|) and ( prj| ) and is the 
relativistic analog of the r-mode. But it is no longer purely axial. 



polytrope with 1.74 solar masses and a radius equal to 12.37 km. The star was very slowly rotating (the ratio between 
its kinetic energy and its mass energy was about 10 -16 ) and we took for a time unit the inverse of the angular velocity. 
The anelastic approximation with the free surface BC was used and to regularize the solution, we added a degenerate 
viscosity with E s — 5. 10~ 5 once the final equations were written in a Newtonian way. We insist on the fact that 
this viscous term does not come from relativistic calculations and just aims at regularizing the solution. In figures 
(f23|), (|24|), (|25|), (26) and (27) are plotted exactly the same quantities as in the Section 0. The conclusions for this 
preliminary relativistic calculation are the same as in the Newtonian case: a polar counter part of the velocity appears 
from the beginning of the evolution and most of the time, the velocity is concentrate near the equator and the surface. 
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FIG. 25: Time evolution of the $ component of velocity in the same calculation as in fig.(|24|). The curve and the associated 
power spectrum show that there is one mode that also appears in the polar part with a frequency very close of the frequency 
of the linear r-mode. See also the fig.(E^). 





FIG. 26: Time evolution of one of the two independent components of the Sij[t] tensor that appears in the RR force. This 
calculation was done during the same run as the results in fig.(|24|) and (|2^). We can see the almost monochromatic associated 
spectrum with exactly the same frequency as the mode that appears in the previous figures. 



VII. CONCLUSION 

This article deals with a new hydrodynamical code based on spectral methods in spherical coordinates. This first 
version of the code was written to study inertial modes in slowly rotating NS, but it can easily be modified to include 
fast rotation. As it was explained in the introduction, the slow rotation limit is a very good approximation for a first 
step. The present version of the code can be used to study both the Newtonian case and the relativistic case with 
the so-called strong Cowling approximation. Here, we have given the main algorithms it uses and shown some tests 
of the linear version. The way we overcame the numerical instability is also explained. 



Furthermore, in order to work with very different time scales, we have introduced an approximation to deal with 
mass (or baryonic number) conservation that proved itself quite robust and useful to study inertial modes: the 
anelastic approximation. Even if this approximation is not a necessary one and can be abandoned in further studies, 
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FIG. 27: # component of the velocity versus the radius for several values of i?. We can see the kind of concentration 
of the motion near the surface and the equatorial plane. This calculation was done with the strong Cowling and anelastic 
approximations and the free surface BC. The background star is a 7 = 2 relativistic rigidly rotating polytrope with 1.74 solar 
masses and a radius equal to 12.37 km. The initial data are the Newtonian linear / = m — 2 r-mode. This figure corresponds 
to the velocity after a time equal to 15 oscillations of the linear I — m — 2 r-mode. 



it is very useful to have an idea of the main properties of the inertial modes in slowly rotating stars. The results 
presented in this article were obtained in the linear case using very basic models of NS, such as the divergence free 
case with a rigid crust (V r | r =i= 0) or a 7 = 2 polytrope with anelastic approximation and free surface. Deeper 
studies in the GR framework with more quantitative results on the effects of EOS, of stratification of the star and 
of BC will follow. Despite the naive aspect of our studies, some common features appeared. Indeed, whatever the 
approximation (divergence free or anelastic) with different BC (rigid crust for the divergence free case and free surface 
for the anelastic approximation), we saw the appearance of a polar part in the mode, as soon as the background is 
Newtonian and differentially rotating, or relativistic and rigidly rotating (in the latter case, this phenomenon is due 
to the frame dragging) . Whether in the case of noisy initial data with a radiation reaction force acting on them, or 
in the case of the free evolution of an initial linear r-mode, from the very beginning the energy of the polar part is 
at least one per cent of the energy of the axial part. This result is compatible with the analytical work by Lockitch 
et al. Q who proved, in the relativistic framework, the existence of a non-vanishing polar part of the velocity in a 
sequence of NS with a barotropic equation of state and decreasing rotation rate. Furthermore, another interesting 
result that we obtained is the "concentration of the motion" near the surface that quickly appears: after less than 15 
periods of the linear r-mode in the calculations done up to now. Adding viscosity has shown that this is a robust and 
physical feature. Moreover, it seems to neither depend on the physical conditions or on the BC. This phenomenon is 
very similar to the results found by Karino et al. in 2001 Q with a eigenvectors calculation. But, if this is really a 
singularity of the derivative of the velocity near the surface, as it seems to be in our time evolutions, it would mean 
that a mode calculation with finite difference schemes can only give quantitative results completely depending on the 
number of points. Indeed, with this scheme, there is an intrinsic numerical viscosity that depends on the resolution. 
Moreover, if this phenomenon is verified in following studies in GR, the stratification of the NS or the physics of the 
crust could be very important to know if inertial modes are relevant for GW emission. 
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The numerical algorithm adopted to solve the NSE (or the EE) is based on the pseudo-spectral methods (PSM) 
widely used in hydro or MHD problems. Before explaining this algorithm in more details, we will begin with a short 
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summary of the PSM in order to make more evident the peculiarity of the solving of vectorial equations like NSE. 
Then, follows a second appendix that deals explicitly with EE and aims at explaining the way we implemented the 
approximation done on the mass conservation equation. 



APPENDIX A: SPIRIT OF THE PSEUDO-SPECTRAL METHODS 

Our group developed algorithms and routines library (see Bonazzola et al. 1990, 1999 (57j],|38]]) 9 allowing us to 
solve partial differential equations (PDE) in different geometries, mainly in domains diffeomorphic to a sphere. First, 
with an example of scalar PDE, we will look at the singularities contained in operators expressed in spherical like 
coordinates [see for instance the different components of the NSE: Eq.@] and at our choice of spectral basis. Then, 
we shall discuss the further difficulties that arise in vectorial PDE and explain the way we overcome them. 



1. The scalar heat equation 



Consider the heat equation: 

dF 

— = aAF + S (Al) 

where AF = ^£ + f §7 + 7^(§^F + fhvf W + sinV f^") ^ s ^ ne Laplacian in spherical coordinates, a the heat 
conductivity supposed to be constant, and S a source term ( tha t may include non-linear terms if they exist). The 
idea of spectral methods is to look for the solution of the Eq.(Al) on the form 

oc 

F[r,4,<p,t]= F nlm [t]H l n [r]KrmL m [ip} (A2) 

n,l,m=Q 

with (0 < r < 1; < <& < w; < <p < 2tt) and where H l n [r] , [&] , L m [<p] are a well chosen complete set of functions. 
The problem is then to find the time evolution of the coefficients F n [ m [t] . 

In a spherical geometry, it is quite natural to choose L m [ip] — exp lmv and K™[d] = where -P™[$] are the 

Legendre functions. With these choices, the Eq.(Al) can be written 

« - a + 2 - l -^F lm [r, t]) = S lm [r, t] (A3) 

where Si m [r, t] are the Fourier-Legendre coefficients of the function S at the radius r and the instant t. 

In order to handle the singularity at r — 0, we shall consider separately the cases 1 = 0,1 = 1 and I > 1: 

- For I = 0, we use the fact that a C 1 function symmetric with respect to the inversion r — > — r has its first 
derivative that vanishes at least as r at r = 0. Therefore, for such a function, the term + is regular. 
Even Chebyshev polynomials T 2n [r] have this property, and the choice H® [r] = T 2n [r] (n > 0) then satisfies the 
regularity conditions. 

- For I = 1, it is almost the same, but the final choice is H^[r] = T2„+i[r]. 

- The case / > 1 is more delicate to handle. Indeed, it can be shown (See Bonazzola et al. 1990 [[57]]) that for 
F[r, -d, (p, t] to be a C°° function, the coefficients Fi m must vanish as r l at r = 0. It means that the Fi m [r, t] are 
symmetric with respect to the inversion r — > — r if I is even, and anti-symmetric in the opposite case. We shall 



In 1980, one of us (S.B.) started to build a library of routines based on spectral methods to solve PDE in different geometries. Today, 
this library contains more than 700 routines written in FORTRAN 70 and 90 languages. These routines are highly hierarchised and 
allow us to assemble codes in modular way. We call this library "Spectra". A part of this library (the highest in the hierarchy) was 
written in C ++ language by J. A. Marck and E. Gourgoulhon in order to allow the use of an object oriented language. This library is 
called "Lorene" . The code described in this section uses the "Spectra" library and is written in FORTRAN 90 language. 
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then distinguish the two cases I even and I odd. 
Case I > 1 even 

The functions H l n [r] = T 2n [r] - T 2n+2 [r] are even and vanish as r at the origin. Therefore the quantity 
(373- + 737 — ii^li )H l n [r] is regular at the origin and is retained. 

Case I > 1 odd 

The functions H l n [r] — (2n + l)T2„+i[r] + (2n — l)T2 n +3[r] are anti-symmetric with respect to the inversion 
r — » — r and vanish as r 3 at r = 0. Therefore the quantity (^3 + 737 — ^^ 1 " > )H l n [r\ is also regular and 
completes our basis. 

Note that by this way, we expand the solution with a set of functions that satisfy minimal conditions of regularity, 
excepted for I = and 1 = 1. It means that they vanish as r 2 or r 3 , in order to make regular any term in the 
equation, instead of vanishing as r l as they should to form a C°° function. Until now, in all the problems that our 



group treated, these minimal regularity conditions were sufficient (see Bonazzola et al. 1999 |38 for more details). 
But we shall see later that they make the Euler equations unstable. 

To end with the scalar case, we shall show how the Eq. ( |A3| ) can be written when the time discretization is performed 
and a 2 nd order implicit scheme used 10 : 

pi+l _ At v^°° A 1 J?i +1 — TP J -^AtfnlX^ 00 A 1 TP 3 4- qi+ 1 / 2 \ f A /I s ! 

r nlm u 2 Z^p=0 ^ L pn r plm ~ r nlm ' y 2 Z^p=0 V plm °nlm J • 

Here A l j n = (f/j- O l H l n ) is the matrix of the operator O l = + 737 — ^75^ • The index j + 1/2 means a quantity 
S computed at the time P + At/2, obtained by extrapolating this quantity using its (known) values at time t 3 and 
P: 5 J+1 / 2 = (3S J — 5 fJ_1 ) /2. The left hand side operator can be easily reduced to a pentadiagonal one, and then 
inverted with a number of arithmetic operations that is proportional to the maximum value N max of N . 

In the case of a space dependent heat conductivity a[r, t9, <p], we can use a semi- implicit method (see Gottlieb et al. 
1977|Q). It consists in solving the equation 



Ktn\ ~ + ^ 2 )f E7=0 KnF^n = (A5) 
Kim + f (« + Er=0 4»^m + At (^™ 2 + (« - O - 6r 2 ) (E^O 4„^ m ) J+1/2 

where a and b are two coefficients that satisfy the condition a + br 2 > a everywhere and that we choose in such a 
way that a + br 2 takes the same values as a at the boundaries. 

The new matrix in the L.H.S. is again a penta diag onal matrix. Boundary conditions are imposed by adding a 
homogeneous solution of the Eq.(A4) or of the Eq.(A5) n . The reader can find in the quoted literature more details 



on the spectral methods, and on how to implement boundary conditions. 

To conclude, note that in the present example, we expand the solution in spherical harmonics. But it turns out 
that in general it is more convenient to use linear combinations of Chebyshev polynomials to treat the expansion 
in 1?. The philosophy of the spectral methods then consists in performing simple operations such as computing 
derivatives, primitives, integrals, multiplications or divisions by r, r 2 , sin$ or cos$ in the coefficients space, and by 
making multiplications of functions in the configuration space. If necessary, the passage to a Legendre representation 
is performed with a matrix multiplication. 



Due to the use of Chebyshev polynomials, the Courant condition s fo r the explicit problem formulation are quite severe. In the present 
case At m ax is proportional to 1/N^ lax (see Gottlieb et al. 1977 |58[), therefore implicit or semi-implicit formulation of the problem is 
required. 

The case a[l, "d, <p] = is called degenerate. In this case no BC are allowed, and a + b must be = 0. 
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2. Analysis in vectorial case 

a. Stability of numerical vectorial PDE 



We have seen how to handle coordinates singularities appearing in scalar PDE when spherical coordinates are 
used. But when treating vectorial PDE, as EE or NSE, the singularities are more malicious: a single look at the 
components of NSE [Eq.(^)] should be sufficient to convince the reader. The singular terms cannot be handled only 
by laying down analytical properties, as it was done in the above scalar equation. Indeed, the dependence existing 
among the different spherical components must be taken into account in order to make singular terms compensate 
each others. A simple example will illustrate this crucial point. 

Consider the following constant divergence free vector V of Cartesian components: V x = 0, V v = 0, V z = 1. Its 
spherical components are V r = - costf; V& = sini?; V v = 0. Then, we have divl? = ^ + *V r + + iff Kj) = 0. 

A small error (for example a round-off error) in the computation of components will obviously forbid an exact 
compensation between the two singular terms 2V r /r and ^(d$ V$ + f^fyK?)- The consequence is the creation of high 
order coefficients and possible instabilities appearing in the iterative process. 



Indeed, we found that the hydro-code for solving linearized EE was unstable. As expected, high frequency 
terms were exploding after few hundred time-steps (few periods of the r-mode) either in the case of the anelastic 
approximation or in the case of the incompressible approximation (divW^ = 0). As explained above, the main reason 
of this instability was the non exact compensation of the singular terms contained in the different source terms in 



the Poisson equations (B7) and (JB9|) (see below). 



To overcome this problem, we define two angular potentials Th and P Q in a such a way that 

W# = d^P a - —„d v T h (A6) 

W v = —8 V P + do T h . 
sm» 

If r P and Th are any arbitrary set of regular functions, we are now sure that the corresponding compo- 
nents W® and W v will make t he d ivergent terms appearing in the vectorial PDE compensate each other. More- 



over, the system of equations (|A6|) can be easily inverted. First, taking the angular divergence of W$ and Wip: 
(^YVi) v W={d§ + cot i?)W# + ^$d v W v gives A# V P = c\yv^ v W where A^=3| + cotddi) + ^z-gd^. If div^VK is 
then expanded in spherical harmonics, we immediately have 

Po,hn = - j^-yj idiv^W^m- (A7) 

It is the same to compute Th- Taking the angular curl of W, we obtain 

A^Th = curl^W (A8) 

where curl^^W= - -^lid^Wv - ^(sintfW,,)). 

To complete this procedure, we have to guarantee that W r can also compensate the singular terms generated 
by the components W$ and W v . A satisfactory way to proceed is to take as a new variable the quantity divVT. 
Indeed, once the quantity g[r, (p] = divVT is known, it turns to be easy to find W r . Bearing in mind that 
divVy = d r W r + 2 W r /r + 1/r 2 div^W", we have, after an expansion in spherical harmonics 12 

1 f 

W rt i m = — u (u ff [u,i?,^], m + l{l + l)P 0iZm [u,tf,v?])du. (A9) 
r Jo 



It is important to take into account the following properties of the spherical harmonics decomposition of the functions P and : for a 
given I the product by r of the component W r i m of W , a regular vectorial function, vanishes as r l ~ 1 . The corresponding poloidal part 
Po,lm behaves in the same way. The associated toroidal part is T„ Ji—i) m and is therefore a regular function. 
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Thus, at each time step, the strategy consists in calculating the potentials P and T/j, and then in ca lcul ating back 



the components W# and W v . The component W r is also computed at each time step by using the Eq.( A9) 



Finally, note that with this procedure, the potentials P Ql Th and the component W r have the correct analytical 
behaviour on the axis e z , i.e they vanish as sin m Thus, the instability of the EE hydro-code was drastically reduced, 
but not completely suppressed. It still kept on growing, but about ten times slower than before. The reason of this 
residual instability was that, at each time step, the round-off errors did not have the correct analytical properties at 
r = 0. Indeed, they did not vanish as r l . Consequently high spatial frequency terms were still generated that did 
not have the good analytical properties and a runaway toward the instability was once again generated. To avoid 
this phenomenon, we project, at each time step, all the scalar quantities (divW^); m [r], P ,i m [r] and T^,;„, [r] on a 
Legendre space P„[r] in such a way to satisfy the analytical conditions. Note that the necessity of this procedure is 
due to the fact that in the EE, solved with spectral methods, no dissipative term (numerical or physical) is present. 
Consequently the numerical instability are not dumped. 

b. Solving the vectorial heat equation 

Now, more details concerning the algorithm to solve vectorial PDE with spectral methods. In solving NSE, we need 
to solve a vectorial heat equation that can be reduced to an equation of the type 

— + V A (//V AB) = J (A10) 
at 

where J is the divergence free term of the source. We remember that it is obtained by the decomposition of the 
source term J in a potential part V</> and a divergence free part J: J = J 4- V</>. 



The Eq.(AlO) can be written in the following form: 



— = /iAB+ (V AB) A Vfi + J. (All) 

From a numerical point of view, the term S — (V A B) A V/i is considered as a source and is computed by a second 
order scheme using its values at the times tj—\ and ij-2- 



Therefore, we have to solve the equation 



no 

— =LiAB + j + S (A12) 



with the condition V • B = 0. 



The technique generally used was already described in Bonazzola et al. 1999 |3q| . Nevertheless, we adopted here a 
slightly different approach which gives more accurate results. 



To solv e the Eq. (A12) it is convenient to introduce the poloidal and toroidal potentials P a and as defined in 



appendix A 2 a . Here we will only consider the case of /x = 1, the semi- implicit scheme which must be used in the 



more general case being straightforward. With a second order scheme, we have in the harmonic representation 

o,lm 2^M dr^ 

pj . \i ( I ( ,> r . . 1 p> ±n ' ) it 

MMm ^ lAL \ 2 \ dr 2 ^ r dr r 2 r o.lm r 2D rlra I +11; 



pj-l-i 1 A + " - o.l-m. , 2 o.lm H'-t-iJ pj-t-i z Pi \ — aiQ^ 

r o ,lm ~ 2^ r I ~d~F 1 h 7~ 37 r^^ojm _ 7* D n m \ — 

' X ( , 2 1(1+1) pj 2 ,5 M . TTJ + 1 /2 



where II is the poloidal component of J + S. Note the presence of a singular term ^B\ +1 in this equation. Bearing in 

DJ+l 



mind that div_B = and that P^\ m vanishes as r l 1 , it is obvious that there should be a compensation between these 
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two terms. Nevertheless, we shall slightly modify the Eq.(Al3) in order to obtain more easily this compensation: 



pj+l 1 At ( d2P ti™ | 2 dP a,lm 1(1-1) pj+l \ _ /A-Ml 



M),im T iit I 2 I dr 2 T r dr r 2 r o,!m I r 2 I "^Jm ^ 1 M),im ^ ^rim T L r 0Um J ^ LL hn 

Once P^ m known, B 3 r ~\ m is obtained as explained above [see Eq.(A9), with g = 0]. 

Finally, the equation for the toroidal part of B reduces to an ordinary scalar heat equation: 
r i+i 1 A / d2r S-i)m 2 dr^(,_ 1)m l(Z-l) - +1 \ 

J /i,(i-l)m "T 2 1 dr 2 r dr r 2 1 h,(l-l)m J ~ \ AL0 > 

rpj I At I - I d2T ^' m _L 2 dr M"» _ ^(^~ 1 ) T j ] , , 9 J + l/2 \ 

J ft,(z-i)m "T" ^ l U dr 2 r dr r 2 h ' lm / C-i)™ J 

where i? is the toroidal component of J + 5. 

APPENDIX B: SOLVING EULER OR NAVIER-STOKES EQUATIONS 

In this section, we shall show how to implement the different schemes proposed in the Section III to handle the 



different time scales appearing in the EE for a slowly rotating star. In all cases, for simplicity, we will only consider 
the linearized equations. First, we will begin with the case where the exact system of equations is to solve. Then we 
shall consider the implementation of physical approximations and finally the problem of BC. 

1. Exact Euler equations 

As usual, we start by doing a decomposition of the velocity vector W in a divergence free component W and a 
potential one: 

W = W + V$, div# = 0. (Bl) 
In a second order time scheme, the divergence free component of the EE is then written as 

W +1 = W + At [-2(1 AW + F]°+p /2 (B2) 

while the potential part is 

&+ 1 = & - At [ip Q + h + <p] j+1/2 . (B3) 



In the Eq.(BS), the subscript DF means the divergence free component of the external force. In the Eq. ( |B3[ ) , 
the subscript O is for the potential component of the Coriolis force 2ft A W, <p is the potential component of some 
other exterior force and h is the variation of the enthalpy. Moreover, we also have the baryonic number (or mass) 
conservation equation: 

h^ = V-At((w r + § ) d -^ + rH Q A^ +1/2 (B4) 

where Hq is the enthalpy of the non perturbed configuration. In this explicit second order scheme, the required 
typical value of the time step At, which guarantees the numerical stability, is determined by the term d r Ho that is 
much larger than the Coriolis term for a slowly rotating NS. 
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The implicit version of the system of equations (B3), ( JB4| ) is the following: 

^j+i = 9 j _ At (i + M) + bo + + </?P' +1/2 ) 
/ji+i = ftj _ i At ((d r W +1 + d r ^ j + Wj +1 + wA d r H + TH (A^ +1 + AW] 



(B5) 



where W} +1 is obtained from the Eq.(B2). Solving the above system of equations reduces to invert, for each value of 
I and m an enneadiagonal (2 N r x 9) matrix (9 diagonals) where N r is the number of degrees of freedom in r. The 
generalization to the NSE is quite obvious and we shall not give more details about it. 



2. Implementation of physical approximations 



As it was already explained in the Section III , for solving the EE or NSE, we chose to use approximations in order 
to better control the results. 



a. The divergence free approximation 



The spirit of this approximation consists in replacing the mass conservation equation (B4) by the condition 

divV^ = (B6) 



or equivalently 5* = 0. The problem is then to find the enthalpy h in a such a way that the condition (B6) is satisfied. 
This can be done easily by solving the Poisson equation 



Ah - divF J+1/2 = 



(B7) 



reached by taking the divergence of the EE. As the gradient of an harm onic function is a divergence free vector, such 
a function can then be added to a particular solution of the Eq.(B7) in order to satisfy the boundary conditions. 

Once h is obtained, the EE has to be solved with F = F — V/i. If viscous terms are present, NSE case, a vectorial 
type equation must be solved. 



Anelastic approximation 



As it was already explained, the anelastic approximation consists in neglecting the term dth in the Eq.(B4). Once 
again, the game consists in finding the enthalpy h in a such a way that the following equation is satisfied: 



TH div W + W r d r H = 0. 



(B8) 



To find h, it is sufficient to derive with respect to time the Eq.(BS), and to replace <9tdivW by its value reached 
from the EE. It gives 



r H Ah + d r H d r h = r Hq divF j+1/2 + Fj +1/2 d r H 
where F contains all the force terms (Coriolis force included). 



(B9) 



The solution of the above equation is achieved with a semi-implicit scheme very similar to the one used to solve 
the Eq.(Al). The problem is then reduced to solve the following Equation (for simplicity, we shall consider only the 
case of a polytropic EOS): 



(7 - l)(a + br 2 )Ah + 2brd r h = H a Ah + d r H d r h + divF j+1/2 + <9 r H F J r +1/2 



(BIO) 



where a and b are t wo co nstants defined as was done in solving the Eq.(Al). The L.H.S. operator can be easily 
inverted and the Eq.( B10| ) can be solved by iteration (see Gourgoulhon et al. 2001 p9[). The iteration at each time 
step being quite time consuming, a more convenient strategy (although less accurate) consists in replacing h in the 
R.H.S. of the Eq.plCl) by h j+1 / 2 . 
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3. Boundary conditions 



It was already explained how to impose BC for divergence free case. For the system of equations (B5) or the 
Eq.(B9), it is worth distinguishing two cases: either the enthalpy vanishes or does not vanish at the boundary of the 
integration domain. Here we will only consider the case of EE. Indeed, for NSE we chose to take a viscosity that 
vanishes at the surface to avoid the need of further BC. 



If Hq \ r —i > 0, the solution is quite easy 
solution that can be used to satisfy one BC. 



the Eq.(B9) or the system of equations (B5) admit an homogeneous 



If ffo(l) vanishes, the system of equations (B5) or the Eq (B9) are degenerate and therefore no BC can be imposed. 
But we shall show that, in this case, the correct BC are automatically satisfied. 



Indeed, consider first the case of the linearized exact system of equations (B5). On the surface of the non perturbed 
star, we have dth + W ■ VHq that is the the correct BC. The surface Hq + h = then defines the profile of the 
perturbed star. To show that the solution in the anelastic approximation also satisfies the BC, it is more convenient 
to first examine the non- linear case. Here, h is not infinitesimal, and the Eq.(B8) must be replaced with the equation 



T(H + h) div W + W ■ V(H + h) = 



that is satisfied if h is a solution of the non linearized Eq (B9): 

T {H +h)Ah- dWF j+1/2 + (V/i - F j+1/2 ) ■ V{H + h j) = 0. 



(Bll) 



(B12) 



Once again, the surface of the star is defined by Hq + h = 0. But from the Eq.( |Bll ) we have the correct BC: 
W ■ V(H + h) \n +h=o = 0. T he su rface of the star is then an unknown quantity that is determined by the relation 
Hq + h = in solving the Eq. ( B12 ) by iteration. This technique was developed in a different context in the already 
quoted paper [Egl. 



In the linear case, we have W r = at r = 1 [cf. Eq.(B8)], and the enthalpy does not vanish on the non-perturbed 
surface of the star. But once again the free surface of the star can be obtained by looking for the position where 
the total enthalpy Hq + h = vanishes. This surface coincides within first order quantities with the non perturbed 
surface r = 1. Note that if 7 > the pressure P\ r= \ oc h 1+1 \ r= i vanishes within terms o(h). 
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